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Nonclassical light from few emitters in a cavity 


D. Pagel, A. Alvermann, and H. Fehske 

Institut fiir Physik, Ernst-Moritz-Arndt-Univer sit at, 17487 Greifswald, Germany 

We study the characteristics of the light generated by few emitters in a cavity at strong light- 
matter coupling. By means of the Glauber ^^^^-function we can identify clearly distinguished pa¬ 
rameter regimes with super-Poissonian and sub-Poissonian photon statistics. We establish a relation 
between the emission characteristics for one and multiple emitters, and explain its origin in terms 
of the photon-dressed emitter states. Cooperative effects lead to the generation of nonclassical light 
already at reduced light-matter coupling if the number of emitters is increased. Our results are 
obtained with a full input-output formalism and master equation valid also at strong light-matter 
coupling. We compare the behavior obtained with and without counter-rotating light-matter inter¬ 
action terms in the Hamiltonian, and find that the generation of nonclassical light is robust against 
such modifications. Finally, we contrast our findings with the predictions of the quantum optical 
master equation and find that it fails entirely at predicting regimes with different photon statistics. 

PACS numbers: 42.50.-p, 42.50.Ar, 42.50.Pq, 03.65.Yz 


I. INTRODUCTION 

Two-level emitters interacting with a cavity photon 
mode are widely studied in quantum optics with respect 
to spontaneous emission and superradiance urn , coop¬ 
er at ivity and lasing m\, as well as the emission of non¬ 
classical light [HHin]- For sufficiently weak light-matter 
coupling, when the photon-dressing of emitter states is 
negligibly small, the emitter-cavity system can be stud¬ 
ied with the quantum optical master equation, usually in 
combination with the rotating-wave approximation m- 
The quantum optical master equation fails at strong 
light-matter coupling, to the extent that it predicts un¬ 
physical emission at zero temperature if the number of 
photons in the ground state is finite [IMl]. 

The correct theoretical description of systems with 
(ultra-)strong light-matter coupling p^bUTS] has attracted 
increasing interest recently [T9]-[24]. Essentially, the 
quantum optical master equation has to be replaced 
by a master equation expressed in the photon-dressed 
emitter eigenstates [24H29]. While the master equation 
remains Markovian, which is justified because of the 
weak emitter-environment and cavity-environment cou¬ 
plings mm, it now requires full diagonalization of the 
interacting emitter-cavity Hamiltonian. Such an equa¬ 
tion was used in recent studies of photon blockade ef¬ 
fects m, spontaneous conversion of virtual to real pho¬ 
tons Eniin], and the emission of nonclassical light from 
a single emitter [22]. 

In this paper we study the emission of few emitters in 
a cavity, with particular focus on the photon statistics 
of the emitted light. Our goal is the characterization 
of temperature and coupling regimes where nonclassi¬ 
cal light [30l [31] is generated. A major result will be 
the identification of two clearly distinguished neighbor¬ 
ing regimes with pronounced sub-Poissonian and super- 
Poissonian photon statistics at strong coupling. 

Our results are obtained with the full input-output for¬ 
malism [ 32 II 35 ] and master equation [24ll29] without fur¬ 
ther approximations. To understand the relevance of the 


different approximations involved in traditional quantum 
optics treatments we make two comparisons. First, we 
compare the results that are obtained when the counter¬ 
rotating light-matter interaction terms are included in 
the Hamiltonian to those when they are dropped. Sec¬ 
ond, we contrast the results obtained with the full master 
equation with results from the quantum optical master 
equation. The latter comparison will clearly show the 
necessity of using the correct master equation already at 
weak coupling if the photon statistics is of interest. This 
issue has been studied conclusively for a single emitter in 
Ref. [ 22 ], which also contains Glauber function plots for 
few emitters in the supplemental material but omits the 
further analysis of the situation that we give here. 


The paper is organized as follows. In Sec.[ll|we intro¬ 
duce the physical situation under study together with the 
master equation used for its analysis. In Sec. |III| we dis¬ 
cuss the emission spectra in relation to the energy spectra 
of the emitter-cavity Hamiltonian, while the statistics of 
the emitted photons is studied in Sec. IV We conclude 
in Sec. |V| The appendices collect further information on 
the theoretical approach. App. ^ gives details of the 
input-output formalism. In App. [^we derive the master 
equation, and give a few analytical results for the photon 
statistics in App. 0 


II. THE PHYSICAL SITUATION 


The interaction of N two-level emitters with a single 
cavity photon mode is described by the Dicke model [36] 

N N 

H = ujcO^a + ujx + g + acr^^) 

i=i i=i 

N 

+ ^ 0 ^ 77 ), ( 1 ) 

i=i 
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where the operator annihilates (creates) a cavity pho¬ 
ton with frequency ujc and ^ is the correspond¬ 

ing lowering (raising) operator for the jth emitter with 
transition energy Ux- Throughout this work, we consider 
the resonant case ujq = uJc = oox- We allow for different 
emitter-photon coupling strengths for the co-rotating (g) 
and counter-rotating {g') interaction terms. Changing g' 
relative to g interpolates between the Tavis-Cummings 
(TC) limit {g' = 0) without and the Dicke limit {g' = g) 
with counter-rotating terms. Both situations can be real¬ 
ized experimentally [13 EH] . The rotating-wave approxi¬ 
mation consists in replacing the Dicke by the TC limit. 

Dissipation arises from the coupling of the emitters and 
the cavity to the environment. For a bosonic environment 
the coupling terms are of the form 

Hi = -\SY,^AK-bl), (2) 

V 

where S' is a (Hermitian) emitter or cavity operator 
and the are bosonic operators for the environment 
photons (at frequencies Uy with coupling constants A^y). 
As the operator S we choose the field operator X = 
—iXo(a — o)) for the coupling of the cavity and the tran¬ 
sition operator cjy^ = i(cr^^ — crff^) for the coupling of 
the jth emitter to the environment. 

At sufficiently weak coupling to the environment, the 
emitter-cavity system density matrix p obeys a Marko¬ 
vian master equation [241429] 

= -i[H, p{t)] - ii e(^) [sls^, Pit)] 

UJ 

UJ 

(3) 

where 

5'u = E (4) 

m,n 

is the projection of S onto transitions between eigenstates 
|m), IjA of with energy difference uunm = En —Em (see 
App. [Bjfor a derivation). For the sake of notational sim¬ 
plicity we state the master equation for a single coupling 
term Q . Multiple coupling terms lead to additional con¬ 
tributions of the same form. 

The functions x(co’) and ^(cj) in Eq. ^ follow from the 
environment spectral function 

Jiu) = 2 tt'^XIS{u! ( 5 ) 


and 

fRer(w + iO+)[n(a;,T) + l] ifa;>0, 

—Rer(—cj + i0+)n(—cj, T) if cj < 0 , 

with the Bose-Einstein distribution function 

Note that in the zero temperature limit n{uj^T) —> 0 
such that the master equation ^ contains only dissipa¬ 
tive terms for transitions |n) -^ym) with positive energy 
^nm > 0 , he., dissipation correctly leads to energy de¬ 
crease. In particular, the problem of unphysical emission 
from the ground state encountered for the quantum op¬ 
tical master equations is resolved. 

In the present work we assume an Ohmic spectral func¬ 
tion 7 c(co’) = yco’/co’o for the cavity-environment coupling, 
and use 7 = in all numerical computations. To 

reduce the number of free parameters we assume the 
same spectral function ^x\uj) = 7 c(<^) for the emitter- 
environment couplings. The respective environment tem¬ 
peratures are also identical. 


A. Solution of the master equation 


As we show in App. the master equation ^ splits 
into two equations of motion 

^Pn,n(^) — ^ ^ X{^kn)En,kPk,k{t) 

kj^n 

^ ^ X{^nk)Ek,nPn,n{t) , (9) 

kj^n 

^Pm,n(^) — i^Em T Pm,n{t) i {jkl 7 ^ Tl) (10) 

for the matrix elements pm,n{t) = (m'|p(t)|n) of the den¬ 
sity operator. In these equations, Sn,k = \{'^\E\k)\‘^ and 


Efi — ^ ^ ^ \x{^nk) T iC(^n/c)] ^/c,n T IT'n . (ff) 


k^n 


The general solution of Eq. (10) is 

Pm,nit) = , (m 7 ^ n) . 


( 12 ) 


Because ReZ^ > 0 for all n, the off-diagonal elements 
of p{t) decay exponentially. Hence, the stationary state 


fulfills 


and its analytical continuation T{uj) into the upper half 
plane, with ^{uj) = Tliiir(±co’ + i0+). Eor a thermal 
environment with inverse temperature (3 = 1/T we get 


xi^) 


7 (w)[n(w,T) + 1] 
7 (-w)n(-w,T) 


if cj > 0 , 
if cj < 0 


(6) 


Pm,n — hm Pm,n{t) — Pn,n^rn,n • (13) 

The diagonal elements p^^ are determined by the sta¬ 
tionary solution of the Pauli master equation If the 
system is coupled to a thermal environment as in Eqs. § 
and 0, the stationary solution of Eq. 0 is the thermal 
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state oc of the system corresponding to the 

temperature T = 1//3 of the environment. 

The emission spectrum and photon statistics can now 
be computed through a standard input-output formal¬ 
ism (see App. 0’ which leads to the projected cavity- 
environment coupling operator 

X_ = -i ^ {En-E^)\m){m\X\n){n\ (14) 

m,n>m 



describing the emission. The correlation functions of X_ 
and characterize the properties of the emit¬ 

ted light. The emission spectrum of the cavity is 

S{u))= lim + r)X_(i))dr, 

t^oo tt Jo 

(15) 

and the second-order Glauber function [39] reads 


( 2 ) ^ {X+{t)X+{t + T)X.{t + T)X.{t)) 

^ ^ {X+{t)X_{t))^ 

Note that evaluation of Eqs. ( [16| ) requires diagonal- 
ization of the Hamiltonian H. 

Because the stationary state from Eq. (13) is di¬ 


agonal in the eigenbasis |n) of we can evaluate the 
r-integration in Eq. ([T^ analytically as 


(17) 


S{u}) = \{m\X-\n)\^p^ 

IT ^' 

m<n 

Iie{Zm + Zn) 


[uj - lm{Zn - Zm)]^ + [Re(Z^ + ZJ]^ ' 

The emission spectrum S(u;) is the sum of Lorentz peaks 
with width Re(Z^+Z^) at the respective transition ener¬ 
gies Im(Z^ — Z^), which according to Eq. (11) are shifted 


relative to the transition energies En — Em of the closed 
system by a Lamb shift that results from coupling to the 
environment. 


B. Quantum optical master equation 

It is instructive to compare the master equation (§ to 
the quantum optical master equation m 

^p(t) = -i[H,pit)] -iYYiSiS^,pit)] (18) 

which is obtained by replacing the projected operators 
Suj with the ‘bare’ operators S± = 

suming x(±Cc;) ^ x±, ^ in the vicinity of a 

typical transition energy uj. Note that 5'+ = —iXoa for 
the cavity-environment coupling and 5'+ = — iaff ^ for the 
emitter-environment coupling. Evidently this approxi¬ 
mation can be valid only for weak light-matter coupling 


FIG. 1. Schematic energy level pattern for the construction 
of the spectrum of the Dicke model. Horizontal arrows depict 
the co-rotating interaction terms in Eq. ^ (coupling constant 
g), diagonal arrows depict the counter-rotating terms (g'). 


g^g' <^uJmUJx, when the dressing of emitter states by cav¬ 
ity photons can be neglected. Because the quantum opti¬ 
cal master equation does not distinguish between energy- 
increasing and energy-decreasing transitions, which are 
equally contained in the unprojected operator S because 
of hermiticity, it can lead to unphysical predictions such 
as emission out of the ground state. Eurthermore, be¬ 
cause failure to observe the above distinction is tanta¬ 
mount to a high-temperature approximation, one will 
expect that the quantum optical master equation fails 
at the prediction of non-thermal photon statistics at low 
temperatures. Therefore, we use the more general master 
equation (§. 


III. THE EMITTED LIGHT 

The first characterization of the light generated in the 
cavity is provided by the emission spectrum. Because 
the emission spectrum depends on the (Lamb-shifted) 
energy spectrum of the Dicke Hamiltonian we start with 
a discussion of the eigenvalues of H for few emitters, 
before we turn to the actual function S{uj) obtained from 
numerical solution of the master equation 


A. Energy spectrum of the Hamiltonian 

To construct the energy spectrum of H we notice first 
that the eigenstates of N uncoupled two-level emitters 
can be classified as angular momentum eigenstates with 
total angular momentum J = A/2, A/2 — 1,... > 0. 
Since H commutes with the total angular momentum op¬ 
erator, states with different J do not mix even at finite 
coupling g^g' ^ 0. Eor fixed J, the Jz quantum num¬ 
ber M can assume the values M = — J, — J -h 1,..., J, 
and a corresponding emitter eigenstate has energy (M -|- 
N/2)uJx‘ Note that for A > 3 the classification in terms 
of J, M is not exhaustive, since different emitter states 
can have identical values. However, these states give the 
same contribution to the emission spectrum. The cavity 
photon eigenstates are Eock states |n) with energies ncjc- 
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FIG. 2. Eigenvalues of H in Eq. 0 for g — g (panels a-c) 
and g — ^ (panels d-f) as functions of the coupling strength 
g. Panels (a,b) show the results for one emitter, whereas the 
number of emitters is A/" = 2 in panels (b,e) and N = 3 in 
panels (c,f). 


For given J we can arrange the eigenstates of the un¬ 
coupled emitter-cavity system as the rungs of a ladder di¬ 
agram as in Fig. Working at resonance ujc = 
the energy (M + N/2 + n)uJo of each state is given by 
the total number of emitter and cavity excitations. The 
co-rotating light-matter interaction terms in H preserve 
the number of excitations and connect states at the same 
energy level (horizontal arrows in Fig. [^. The counter¬ 
rotating terms change the number of excitations by two 
(diagonal arrows in Fig.[^. This simple scheme explains 
many properties of the energy spectra of H shown in Fig. 

El 

For A/' = litisJ=l/2, and we recover the Jaynes- 
Cummings ladder [40] for = 0. The lowest level (M = 
— 1/2, n = 0) does not couple to any other state and 
hence leads to the ^-independent eigenvalue zero of H [see 
Fig-Hd)]. Every other level consists of two ladder rungs. 
They lead to the eigenvalues nuj^ ± ^/ng^ for n > 1. For 
g' = g corrections arise from coupling between states at 
different height in the ladder but the energy level pattern 
remains discernible [Fig. ^a)]. 

For N = 2 [see Fig.^b,e)] we have either triplet 
( J = 1) or singlet ( J = 0) emitter states. For g' = 0, the 
triplet states lead to the eigenvalue zero (n = 0), the two 
eigenvalues ujq ± ^/2g {n = 1), and the three eigenvalues 
ncjo, nujQ ± \/2^j2n — Ig for n > 2. The singlet states do 
not couple with each other and lead to the ^-independent 
eigenvalues (n + l)co’o for n > 0. It follows that the eigen¬ 
values nujo for n > 2 are twofold degenerate (one triplet, 
one singlet state). This degeneracy is lifted for g' = g^ 
but the energies of the singlet states remain fixed. 

For N = 3 we have quadruplet (J = 3/2) and dou¬ 
blet (J = 1/2) emitter states. The ladder scheme for the 
doublet is equal to that for N = 1 and hence leads to 
the same energy spectrum, apart from the fact that all 



0) / COq CO / COq 


FIG. 3. Emission spectra S(uj) for iV = 2 emitters. Panels 
(i-iii) show results for g' — g^ panels (iv-vi) for = 0. The 
emitter-cavity coupling strength and the environment temper¬ 
ature are g = O.Scjq, T = O.OTc^o in panels (i, iv), g = O.Tcjq, 
T = 0.23ct;o in panels (ii,v), and g = O.Scjo, T = 0.1cc;o in 
panels (iii,vi). 


energies are shifted up by ujq when going from N = 1 to 
N = 3. Notice that the doublet states are two-fold de¬ 
generate, because the angular momentum classification of 
the emitter states is not unique in this case. The quadru¬ 
plet states lead to one (starting at zero for ^ = 0), two 
(at cjo), three (at 2co’o), and four (at nujo with n > 3) 
additional eigenvalues in Figs. if). Because of the close 
vicinity of many states in the energy spectrum the correc¬ 
tions resulting from the counter-rotating terms for g' = g 
are large. This trend continues if N is increased further. 


B. The emission spectrum 

In Fig. we show the emission spectrum S{uj) for 
N = 2 emitters at different coupling strength g and en¬ 
vironment temperature T. These data, as well as those 
for the Glauber function g^‘^\t) shown later, have been 
computed with a maximal number of 10^ cavity photons 
in the numerical diagonalization of i7, which is sufficient 
for the given parameter combinations. 

For low temperatures T <C co’o, only the first possible 
transition into the ground state contributes to the emis¬ 
sion spectrum. It leads to the single peak in panels (i) 
and (iii) of Fig. With increasing temperature transi¬ 
tions involving higher excited states begin to contribute. 
For example the two peaks in panels (ii), (iv) correspond 
to the transition from the 2nd to the 1st excited state 
and from the 3rd excited to the ground state. As could 
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FIG. 4. (Color online) Expectation value (X-^X-) for N = 
2 emitters, as a function of temperature T and coupling 
strength g. Panel (a) shows the result for g' — panel (b) 
for g — Crosses mark the parameters used in Fig. 


FIG. 5. (Color online) Glauber function ^^^^(0) at zero time- 
delay for one emitter [N — 1), as a function of temperature 
T and coupling strength g. Panel (a) shows the result for 
9 — 9^ panel (b) for g' — 0. Note that all values ^'^^^(0) > 4 
are assigned the same (dark-red) color in the density plots. 


be deduced already from panels (b), (e) in Fig. the 
transitions tend to have smaller energies in the TC limit 
than in the Dicke limit, which leads to the red shift of the 
emission peaks in panel (v) relative to those in panel (ii). 
However, at not too strong coupling the low-lying states 
still have comparable energies, and the emission spectra 
look similar. The situation changes at ultrastrong cou¬ 
pling when the co-rotating and counter-rotating terms 
are of equal magnitude (panels (hi), (vi)). In addition to 
the markedly different peak energies the peak height has 
now decreased by two orders of magnitude in the Dicke 
limit, but not in the TC limit. 

The decrease of peak height can be recognized in the 
cj-integrated emission spectrum 

shown in Fig. The equality with the given expecta¬ 
tion value follows directly from Eq. 0- Only in the 
Dicke limit, but not in the TC limit, the total emission 
becomes small again at ultrastrong coupling and low tem¬ 
peratures. Still, one sees that both plots agree nicely for 
not too strong coupling {g/uj < 0.5). This observation 
sets the upper limit of the coupling strength (here, for 
N = 2 emitters) below which the presence or absence 
of counter-rotating interaction terms does not affect the 
light emission significantly. We will find the same behav¬ 
ior for the Glauber function. 


IV. NONCLASSICAL LIGHT 

A basic decision on the possible generation of nonclas- 
sical light is possible with the Glauber function ^^^^(0) 
at zero time delay. For ^^^^(0) = 1 the emitted pho¬ 
tons have a Poissonian distribution, while ^'^^^(O) > 1 
indicates super-Poissonian statistics. Thermal light has 
^^^^(0) = 2. By contrast, ^^^^(0) < 1 indicates nonclassi- 
cal light with sub-Poissonian photon statistics. Further 


information on photon (anti-)bunching is provided by the 
full time-dependent function g^‘^\t). 


A. Photon statistics for one emitter 

The Glauber function ^'^^^(O) for one emitter [N = 1) 
is shown in Fig. Two distinct regions can be identi¬ 
fied in the Dicke limit in panel (a) (where g' = g). A 
triangular region with ^^^^(0) < 1, which stretches out 
along the vertical axis, indicates the emission of nonclas- 
sical light with sub-Poissonian photon statistics at low 
temperatures and moderate-to-strong light-matter cou¬ 
pling. It lies below an elongated region with strongly 
super-Poissonian photon statistics (^'^^^(O) ^ 2) at larger 
coupling, which extends diagonally towards higher tem¬ 
peratures. Both regions are embedded in the background 
of thermal light with ^^^^(0) ~ 2. The situation is dis¬ 
tinctly different in the TC limit {g' = 0) in panel (b), 
where the super-Poissonian region is pushed back in favor 
of a second sub-Poissonian region that continues towards 
ultrastrong coupling. Note, however, that the emission 
of nonclassical light in the first sub-Poissonian region is 
observed equally in both limits. 


B. Photon statistics for few emitters 

The distinctive features of the Glauber function persist 
for multiple emitters (see Fig. [^, but the regions are 
shifted to smaller couplings g as the number of emitters 
increases from one to three. 

The obvious similarity between ^^^^(0) for N = 1,2,3 
emitters visible in Figs. [^can be expressed as an ap¬ 
proximate relation between the respective emitter-cavity 
coupling g. In the Dicke limit {g' = g) we find that 
the features of ^^^^(0) are closely reproduced under the 
scaling g oc 1/N. In the TC limit {g' = 0) features 
are reproduced under the scaling g oc 1/^/N. Interest- 
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FIG. 6. (Color online) Glauber function at zero time- 

delay as a function of temperature T and coupling strength 
for iV = 2 emitters (top panels a,c) and N = 3 emitters 
(bottom panels b,d), The left panels (a,b) show results for 
g' = g, whereas = 0 in panels (c,d). Note that all val¬ 
ues ^^^^(0) > 4 are assigned the same (dark-red) color in the 
density plots. 


ingly, the proper scaling of g depends on the presence of 
counter-rotating interaction terms in the Hamiltonian. 
This difference is in contrast to the semiclassical theory 
where the mean cavity photon number in the steady state 
scales oc N both in the Dicke and TC limit. Not surpris¬ 
ingly, the Glauber function ^*^^^(0) is more sensitive to 
the details of light-matter coupling than the semiclassi¬ 
cal theory that neglects quantum correlations in favor of 
a mean-field approximation. 

Our arguments in favor of the above scaling relations 
depend on several observations, which we now develop 
for the TC limit {g' = 0). Without counter-rotating in¬ 
teraction terms the Hamiltonian Ff commutes with the 
operator Nt = a^a + ^f=i cr^^cr^l\ which counts the to¬ 
tal number of excitations. Hence, H is block diagonal 
with blocks of the form ntOJoI-\- gC, where rit denotes the 
eigenvalue of Nt^ I is the identity matrix, and the matrix 
block C contains the ^-independent matrix elements of 
the co-rot at ing interaction terms in H. From this form 
of the blocks it is evident that the eigenvectors of H do 
not depend on i.e., the matrix elements of X± that 
enter Eq. (16) are constant. The dependence of ^^^^(0) 
on g results from the eigenvalues only, which determine 
the occupation of the states in the stationary (thermal) 
state and the prefactors of X±. If we can show that the 
eigenvalues scale approximately as g^/N the above rela¬ 
tion follows. 

Let us focus on the low lying states that give the domi¬ 
nant contribution in the interesting temperature regimes. 


M =-J+2 
M =-J+1 
M = -J 






n = 0 1 2 


FIG. 7. Schematic energy level pattern. Horizontal arrows 
depict the co-rotating interaction terms in Eq. 0 (coupling 
constant g ), diagonal arrows illustrate the action of the cavity 
photon annihilation operator a. 


These states can be found in the ladder diagram of states 
in Eig. They must be connected to the ground state 
at energy zero by a diagonal arrow that gives the action 
of the operator a, i.e., of X_. 


Eor the denominator (X+X_) of ^^^^(0) from Eq. (16) 
states contribute which are separated by one vertical step 
in the ladder diagram. The energy of the most relevant 
first excited state is given by Ei = ujq^ qVN^ which has 
the postulated scaling. This scaling of the first excited 
state for few emitters has been verified experimentally in 
Ref. [41]. 


Eor the numerator of ^^^^(0), where 

each operator appears twice, states contribute which are 
separated by two vertical steps on the ladder. Now the 
second excited state is most relevant, which is the lin¬ 
ear combination of the two {N = 1) or three {N > 2) 
vertical rungs that occur for rif = 2 excitations. The cor¬ 
responding 2 X 2 or 3 X 3 matrix from the above block 
decomposition of H is 


/ 2wo V2g\ 
\V2g 2 loo ) ’ 


/ 2uo \f2g 0 \ 

^^2g 2wo ^g . (20) 

V 0 ^/^g 2wo / 


Diagonalization gives the energies E-i = 2wo ± ^g for 
V = 1, while G {2wo, 2(jJq ± \/2\/N + I 5 } for N >2. 
With the approximation ^/N + 1 « s/N, which is good 
enough for a rule of thumb, this is again the postulated 
scaling. Put together, the energies that enter the com¬ 
putation of 5 '^^^( 0 ) scale roughly as g\/N^ which con¬ 
cludes our argument in favor of the observed relation 
oc l/v/]V” in the TC limit. 

In the Dicke limit g' = g the block decomposition of 
H is not possible because of the counter-rotating interac¬ 
tion terms. The eigenvectors of H now depend on and 
the previous argument cannot be easily translated. How¬ 
ever, inspection of the energy spectra in Eig. [^strongly 
suggests that the observed relation is still related to an 
approximate relation between the eigenvalues of H for 
different N^ now with the scaling g oc 1/N. 
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FIG. 8. (Color online) Glauber function computed 

with the quantum optical master equation ( |18[ ), shown as a 
function of temperature T and coupling strength g for g — g. 
Panel (a) gives the result for iV = 1 emitter, whereas iV = 2 
in panel (b). 


C. Photon statistics from the quantum optical 
master equation 


Results for the Glauber function obtained with the 
quantum optical master equation (18) are shown in Fig.|^ 
in the Dicke limit g' = g. In stark contrast to the results 
from Figs.[^[^the quantum optical master equation does 
not predict the emission of nonclassical light with sub- 
Poissonian photon statistics in any part of the parameter 
space. The situation does not improve in the TC limit 
= 0 where [i7, Nt] = 0 and the quantum optical master 
equation gives the stationary (thermal) state oc 
leading to ^'^^^(O) = {a)a)aa)/{a^a)‘^ = 2 independent of 
the number of emitters N^ the coupling strength g, or 
the temperature T, thus always predicting the emission 
of thermal light. While it may not be surprising that the 
quantum optical master equation fails, because the weak 
coupling condition g uj^.c is not satisfied, it is remark¬ 
able that it fails to capture any features from the previous 
Glauber function plots in Figs. This failure high¬ 

lights the importance of using the correct master equa¬ 
tion not only for strong light-matter coupling but also 
if one is interested in properties following from higher- 
order correlation functions, such as the photon statistics 
obtained from the second order Glauber function. 


D. Photon bunching and antibunching 

A further property to distinguish classical and nonclas¬ 
sical light is the time-coincidence statistics of the emitted 
photons, which can be deduced from the time-dependent 
Glauber function g^‘^\t). For classical light, has 

a non-positive initial slope at t = 0. This indicates pho¬ 
ton bunching, i.e., that the probability of observing two 
photons at equal times is larger than the probability of 
observing them at different times. Conversely, a posi¬ 
tive slope indicates photon antibunching, which is pos¬ 
sible only for nonclassical light. In the long-time limit. 



FIG. 9. Glauber function g^‘^\t) as a function of time for 
N — 2 emitters. Panels (i-iii) show the results for g — g^ 
whereas in panels (iv-vi) g — 0. The emitter-cavity coupling 
strength and the bath temperature are g — 0.5ct;o, T — 0.07cc;o 
in panels (i, iv), ^ = 0.7ct;o, T = 0.23ct;o in panels (ii,v), and 
g — 0.8ct;o, T = O.lcjo in panels (iii,vi). 


limt-s-oo = 1 in all cases. 

In Fig.we plot g^‘^^ (t) for the parameter combinations 
marked in the two upper panels in Fig. We see that 
g^‘^^ (t) is always a strictly monotonic function of t. There¬ 
fore, in the present situation photon bunching and an¬ 
tibunching coincide precisely with super-Poissonian and 
sub-Poissonian photon statistics. Only if 1 <5(2)(0)<2 
in panel (ii) the function oscillates slightly, but 

the overall decay is still indicative of photon bunching. 


V. CONCLUSIONS 


Our analysis of the light generated by few emitters in 
a cavity reveals a non-trivial dependence of the photon 
statistics on the light-matter coupling and temperature. 
Clearly identifiable parameters regimes with sub- and 
super-Poissonian photon statistics appear at strong and 
ultrastrong coupling, and lie immediately next to each 
other. Tuning of the light-matter coupling or change of 
the temperature can thus have a tremendous effect on 
the photon statistics. As a general trend we find strong 
signatures of nonclassical light at strong coupling. Ther¬ 
mal photon statistics, on the other hand, requires weak 
coupling or high temperatures: It is the exception rather 
than the rule at low temperatures. 

The photon statistics, and to a lesser degree also the 
total emission, is strongly influenced by the presence 
of counter-rotating light-matter interaction terms in the 
Hamiltonian. These terms are responsible for the preva- 

























lence of super-Poissonian over sub-Poissonian light at 
ultrastrong coupling. Not surprisingly, the convenient 
rotating-wave approximation (i.e., identihcation of the 
Dicke by the TC limit) gives the wrong prediction when 
the coupling becomes too large. Nevertheless, the sce¬ 
narios with and without counter-rotating terms are sur¬ 
prisingly similar at not too strong coupling, which shows 
that generation of nonclassical light is not a peculiar ef¬ 
fect arising from the hne-tuning of interaction terms in 
the Hamiltonian but a rather robust feature. 

We have provided an approximate rule to relate the 
emission of few emitters to the emission of a single emit¬ 
ter, under appropriate scaling of the coupling constant. 
In accordance with this rule, the features of the Glauber 
function observed for one emitter occur at comparably 
smaller values of the individual emitter-cavity coupling in 
the case of a few emitters. The reason is that ah emitters 
interact with the same cavity mode, which magnihes the 
effects of resonant emission and (re-)absorption of cavity 
photons. Broadly speaking, generation of nonclassical 
light is easier with more emitters because the required 
coupling of each individual emitter to the cavity mode 
can be reduced. 

Our analysis of strong light-matter coupling required 
use of the full input-output formalism and of the full 
master equation, which carefully distinguishes between 
transitions at different energies. If this correct treatment 
is replaced by the standard quantum optical master equa¬ 
tion results change completely. Especially, the prediction 
of nonclassical light does not survive the additional ap¬ 
proximations made in the replacement. While the quan¬ 
tum optical master equation could not be expected to 
work at strong coupling, its outright failure at describ¬ 
ing any of the distinctive features observed in the photon 
statistics shows that use of the right master equation is 
essential in ah situations, perhaps apart from extremely 
weak coupling. The price one has to pay is full diagonal- 
ization of the Hamiltonian. 

We here focus on the system at thermal equilibrium. 
Future work should address emission if the system is 
driven coherently through external photon sources. This 
will require addition of explicitly time-dependent peri¬ 
odic terms to the Hamiltonian, and thus combination 
of the present master equation with the Floquet formal¬ 
ism. By contrast, a perturbative expansion in the driving 
strength is sufficient only for weak off-resonant driving, 
but then the possible new effects would be weak too. 
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Appendix A: The input-output formalism 


We follow standard input-output theory [H sg. 
The interaction Hamiltonian in Eq. © for the cavity- 
environment coupling in the continuum limit is 

Hi = -iX J D{oj)X{oj){b^ - bl)doj, (Al) 

where D{uj) is the environment density of states 
and Ac(cj) the cavity-environment coupling function, 
i.e., the environment spectral function is = 

271D{uj)Xc{ojd. Hi together with the free Hamiltonian 
j D{uj)ujhljh^^duj of the environment photons and the 
commutator — uj') lead to the equation 

of motion 


buj = -iojbuj + X{uj)X (A2) 

for the held quadratures of the environment. For to < 
t < ti, the formal solution of Eq. (A2) is 

b^{t) = + \{co) t 

Jtn 




to 

ti 


- A(w) J . 

(A3) 


We dehne input (output) held operators 

^in(out)(^) = J D{Lo)X{Lo)e~^^^*~*°<^^'>'’buj{to(l))dLO (A4) 

and make use of the spectral function 7c(<^) = jujjujo to 
obtain the input-output relation 

bout (t) = bin (t) + i—X-(t), ( A5) 


where X_ denotes the positive frequency component of 
X, i.e., X_ acts as a lowering operator. The explicit 
dehnition of in the system-energy eigenbasis is given 
in Eq. (14). 


Appendix B: The Markovian master equation 

We consider the dissipative dynamics of the system 
density matrix in the weak system-environment coupling 
limit. For strong coupling within the system the quan¬ 
tum optical master equation predicts unphysical emis¬ 
sion from the ground state El- Going one step back in 
the derivation of the quantum optical master equation, 
the second-order time-convolutionless projection opera¬ 
tor method [25] gives a time-local master equation lead¬ 
ing to consistent results including the counter-rotating 
terms mmi]. Nevertheless, this master equation does 
in general not generate positive dynamics gailSl. This 
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problem was resolved by a recently derived master equa¬ 
tion in the system eigenbasis [24H29] and we here reca¬ 
pitulate its derivation. 

The total Hamiltonian is the sum of the contribution of 
the system, the contribution of the reservoir, and 
the interaction Hj. We note that the interaction Hamil¬ 
tonian in Eq. ([^ is of the general form = SR^ where 
S {R) is a Hermitian system (reservoir) operator. A more 
general coupling Hj = SnRn can also be considered, 
but leads to the same qualitative results. The dynam¬ 
ics of the density operator prit) of the total system in 
the interaction picture is described by the von Neumann 
equation 


= -i[Hi{t),pTit)] . (Bl) 

As a notational convenience, we mark operators in the 
interaction picture with a hat. The interaction Hamilto¬ 
nian and the density operator in the interaction picture 
are defined as 

Prit) = U^{t,0)pT{t)Uo{t,0) , (B2) 

Hiit) = Ul{t,0)HiUoit,0), (B3) 

where the time evolution operator of the uncoupled sys¬ 
tem and reservoir is 


Uo{t, s) = . (b 4) 

In the limit of weak system-reservoir coupling several 
approximations are performed. First of all, within the 
Born approximation initial factorization of the density 
operator is assumed, Pt(0 ) = p{0)pR, and the back- 
action of the system onto the reservoir is neglected, 
pT{t) = p{t)pR. Secondly, the Markov approximation 
is performed by replacing p(r) at retarded times r with 
p{t) at the local time t. In the third place, assuming that 
the reservoir correlation time is small compared to the 
relaxation time of the system, the time integration is ex¬ 
tended to infinity to arrive at the Born-Markov equation 
of motion 

d . . 

= - J T^R{[Hiit), [Hi{t - t), p{t)pR\] }dT, 

(B5) 

where Tri^{-} denotes the partial trace over the reservoir 
degrees of freedom and {R) = 0 is assumed. We further 
assume a thermal reservoir state pR oc and define 

the reservoir correlation function 

C{t) = = C(-r)* (B6) 

to evaluate the traces in Eq. ( |B5[ ). This yields the master 
equation 

d . 

= [S{t-T)p{t),S{t)]C{T)dT + }i.c., (B7) 

where H.c. denotes the Hermitian conjugate. 


We introduce the transition operators in Eq. that 
are the discrete Fourier components of the interaction 
picture 5'(t), i.e., 

= (B8) 


Equivalently, = —ooS^. In addition, we introduce 

the even and odd Fourier transforms of the reservoir cor¬ 
relation function 

/ oo 

C'(r)e-"dr = x(a;)* , (B9) 

-OO 

C(u;) = l r C{t) sgn(r)e-"dT = . (BIO) 

^ J — OO 

For a thermal photon reservoir with spectral function 
7 (co’) the functions x(cj) and ^{uo) are given in Eqs. ^ 
and 0- With these definitions we find 

UJ,UJ' 

+H.C.. (Bll) 


Eq. (Bll) is the standard Born-Markov master equation 


in the system energy-eigenbasis. It contains the dissipa¬ 
tive parts proportional to x(ci;) and the L amb-shift terms 
proportional to ^(co’). Because Eq. (Bll) is not of Lind- 
blad type, it does, in general, not preserve the positivity 
of the density operator. 

Inspecting Eq. (Bll) we recognize that it contains os¬ 
cillating terms proportional to If we assume 

that the relaxation of the system is slow compared with 
all oscillations we can neglect the contribution 

from terms with uj' ^ uj. This approximation is called 
secular or rotating-wave approximation and the master 
equation in the Schrodinger picture simplifies to the re¬ 
sult given in Eq. §• This equation is the Lindblad mas¬ 
ter equation that includes the Lamb shift of the unper¬ 
turbed system energies En as well as reservoir induced 
dissipation effects to lowest order in the system-reservoir 
interaction strength. 

As is already known in the literature, special care has 
to be taken if the spectrum of H is degenerate [UlIlT]. 
But even if the eigenvalues are non-degenerate we 
may have situations where energy differences are degen¬ 
erate, i.e. En — Em = Ek — El ioT n ^ m ^ k ^ 1. The 
consequences of these two different types of degeneracy 
can be understood when we decompose the density ma¬ 
trix into blocks. In particular, we write (Smn) for 
the matrix containing the elements {k\p\l) {{k\^)) with 
Ek = En and Ei = Em- The master equation ^ in this 
block notation reads 


~ Xi^k Em)^mkPkl{^)^, 

k,l 

1 


ln^Ek-Em,Ei-En 


2 ^ 

k 




- EkYp, 


(i)SLs 


ipi^En 




mk 


nk^kn 


^kmPmni^) 5 (B12) 
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where the summations are over different system energies 
only, and the complex function ip{x) = x(x) +i^(x) is in¬ 
troduced. We see that the last two lines in this equation 
are block-diagonal. For m = n, the Kronecker-delta in 
the first line evaluates to such that diagonal blocks 

couple to diagonal blocks, only. For m ^ the first line 
contains terms with k ^ I only, such that non-diagonal 
blocks do not couple to diagonal ones. Nevertheless, a 
non-diagonal block couples to another non-diagonal 
block p^i with k ^ I ^ m ^ n li the respective transi¬ 
tion energies are degenerate. Thus, energy level degener¬ 
acy introduces a block structure implying that a diagonal 
density matrix element couples to non-diagonal elements 
within diagonal blocks whereas energy transition degen¬ 
eracy leads to a coupling of non-diagonal blocks to differ¬ 
ent non-diagonal blocks. We remark that both subtleties 
have their origin in the rotating wave approximation. On 
the one hand, this approximation leads to the Lindblad 
structure of Eq. On the other hand, it results in strict 
Kronecker delta’s between the two transition energies uj' 
and UJ. 


Consider a situation where each degeneracy in the 
spectrum of H as well as in their differences is lifted by a 
small e parameter. Then, each block contains a single el¬ 
ement only, implying that the equations for the diagonal 
density matrix elements no longer couple to non-diagonal 
elements. In addition, any non-diagonal element of the 
density matrix evolves independently from all other ele¬ 
ments. This behavior does not change when we let each 
e ^ 0. In this limit, the equations become independent 
of the e parameters, but are different from the e = 0 case. 
In particular, for every non-zero e ^ 0 we get the two 
equations © and ( [1Q| ) for the diagonal and non-diagonal 
density matrix elements. 

We remark that in real physical systems one will never 
have perfectly equal or equidistant energies because each 
small perturbation will lift the degeneracies. In the theo¬ 
retical description we may argue that the Lamb shift lifts 
degeneracies. Nevertheless, we have to keep in mind, that 
with Eqs. © to ( [IT] ) we can not study effects that rely 
on degenerate energies or degenerate transitions, e.g. the 
perfectly harmonic oscillator or a system composed of 
completely uncoupled identical subsystems are not cor¬ 
rectly described. 


Appendix C: Analytical results in the 
Tavis-Cummings limit 


In this section we derive analytical results for the 
Glauber ^^^^(O)-function in the TC limit {g' = 0) for 
a single emitter {N = 1). 


According to the argumentation in Sec. |IVB] the dom¬ 
inant contribution to the denominator of^^(0) in Eq. 
(16) at low temperatures is that of the first excited state 
with energy Ei = ujq — g. Specifically, the denominator 



FIG. 10. (Color online) Analytical results for the Glauber 
function for one emitter (N = 1), as a function of tem¬ 

perature T and coupling strength g in the TC limit. Panel (a) 
shows the result following from Eqs. (|C1[) and (C2), whereas 


in panel (b) the result is refined for ^ ^ a;o as explained in 
the text. 


(X+A_) is approximated by 

(Cl) 

In this expression the exponential ig the ther¬ 

mal population of the first excited state and the prefactor 
(cjo ~ is the squared transition matrix element of 

X- between the first excited state and the ground state. 

The most relevant state for the numerator of 
at low temperatures is the lowest eigenstate with energy 
E 2 = 2co’o ~ of fho 2x2 matrix given in Eq. ( [^ . 
To evaluate the matrix elements of the operators X± we 
have to consider the four possible transition sequences 
|2,—) —y |l,dz) —y |0) —y |l,dz) —y |2 ,—), where |0) de¬ 
notes the ground state and |n, ±) are the two eigenstates 
with energies E^ = nuJoE^/ng. This yields the expression 

+ i [wo - (V2 - 1)5] {ujI - g^) [wo - (V2 + 1)5] 

+ 5:^ [a;o - (V2 + l)g]\coo + gfy-Pi^-o-V2g) 

(C2) 

approximating the numerator ofg^^^(O). _ 

The results belonging to Eqs. and ( |C2[ ) are plotted 
in Eig.j^a). Compared to the exact numerical results in 
Eig. I^l^a good agreement appears for 0.2 < g/uJo < 0.4 
and low temperatures. Eor high values of g ^ O.dcjo the 
first excited state becomes closer and closer to the ground 
state such that the finite temperature leads to significant 
contributions from transitions not involving the ground 
state. Eor this reason, the upper part of Eig.jl^a) is not 
well reproduced. In contrast, the lower part of Eig.[lQ|^a) 
is not in accordance with the exact numerical results be¬ 
cause our assumption of low temperatures T g < ujq 
does not include the limit ^ > 0. 
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To improve the results in regions with g ujq we ad¬ 
ditionally have to take into account the transition se¬ 
quences | 1 ,+) — 7 ^ | 0 ) ^ II 5 +) for fho denominator and 
|2, T) — y |l,zlz) — y | 0 ) — y |l,dz) — y |2,-t-) for the nu¬ 
merator. The result is shown in Fig. Il^b), where the 
agreement to the exact results in Fig.j^b) is now very 
good for all temperatures and g < O.dcJo. Note that at 


g ^ OAujq a crossing of the eigenvalues of the second and 
third excited state occurs, as can be seen in Fig. |^d). 
This indicates that the role of these states in the cal¬ 
culation of is interchanged. The impact of these 

eigenvalue-crossings would analytically be reproduced if 
we include contributions to from higher excited 

states. 
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